Method for computing spherical conformal and riemann mapping

ABSTRACT

A classical way of finding the harmonic map is to minimize the harmonic energy by the time evolution of the solution of a nonlinear heat diffusion equation. To arrive at the desired harmonic map, which is a steady state of this equation, an efficient quasi-implicit Euler method (QIEM) is revealed and its convergence under some simplifications is analyzed. It is difficult to find the stability region of the time steps if the initial map is not close to the steady state solution. A two-phase approach for the quasi-implicit Euler method (QIEM) is disclosed to overcome this drawback. In order to accelerate the convergence, a variant time step scheme and a heuristic method used to determine an initial time step have been developed. Numerical results clearly demonstrate that the present method far computing the spherical conformal and Riemnann mappings achieves high performance.

BACKGROUND OF THE INVENTION

Field of the Invention

The present invention relates to a computing method, especially to a method for computing spherical conformal and Riemann mappings applied to brain mapping, surface classification and global surface parameterizations.

Descriptions of Related Art

Conformal surface parameterizations have been studied intensively, and most works deal with genus zero surfaces. The basic approaches are harmonic energy minimization, Cauchy-Riemann equation approximation, Laplacian operator linearization, angle-based flattening method and circle packing, among others. In harmonic energy minimization, a discrete harmonic map is introduced to approximate the continuous harmonic map by minimizing a metric dispersion criterion. Due to the conformal nature of harmonic maps from a genus zero closed surface to the unit sphere, Gu and Yau et. al proposed a nonlinear optimization method for genus zero closed surface by minimizing the harmonic energy iteratively on the unit sphere until convergence to a harmonic map. The method has been applied to brain mapping, surface classification and global surface parameterizations.

The evolution of computing a conformal map f from genus zero closed surfaces to the unit sphere is carried out by a nonlinear heat diffusion equation:

$\begin{matrix} {{\frac{df}{dt} = {{- \Delta}\; f}},} & (1) \end{matrix}$

with f to be constrained on the unit sphere. The explicit (forward) Euler method has been used to produce a steady-state solution of the nonlinear heat diffusion equation (1) in some researches. The explicit scheme is attractive for its simplicity. Unfortunately, it is known to have a small stability region that leads to extremely small time steps. While the implicit (backward) Euler method has a much larger stability region, it involves nonlinear systems. In this step, it is crucial to have a successful iterative method and a good way to produce the associated initial guess when the time step is relatively large. To tackle this challenging problem, the semi-implicit Euler methods have been proposed to simplify the nonlinearity of the systems and improve the computational cost in each iteration.

SUMMARY OF THE INVENTION

Therefore it is a primary object of the present invention to provide a method for computing spherical conformal and Riemann mappings that solves a nonlinear heat diffusion equation by a two-phase approach for the underlying quasi-implicit Euler method (QIEM). Phase-I QIEM is used to quickly find an approximate solution which is close to the steady state solution. Then Phase-II QIEM is applied to compute the steady state solution using the approximate solution produced by Phase-I QIEM. The advantages of the two-phase QIEM are that it only needs to solve a linear system and allows a large time step in each iteration. For the iterative methods of solving the steady state ordinary differential equation (ODE) systems, the adaptive methods for controlling the time step in each iteration are proposed to accelerate its convergence. The present invention not only uses a heuristic method to estimate the initial time step but also develops an adaptive method to control the time step so that the two-phase QIEM possesses high performance. Promising numerical results illustrate the efficiency and stability of the present method.

In order to achieve the above objects, a method for computing spherical conformal and Riemann mappings according to the present invention includes the following steps. First carry out evolution of computing a conformal map f from a genus zero closed surface to the unit sphere (the spherical conformal mapping) as well as from a simply connected surface with a single boundary to a 2D disk (the Riemann mapping) by a nonlinear heat diffusion equation. Then solve the nonlinear heat diffusion equation by a quasi-implicit Euler method (QIEM). Next analyze convergence of the QIEM under some simplifications. Lastly accelerate the convergence of the QIEM by using a two-phase approach for the quasi-implicit Euler method to estimate an initial time step and an adaptive method to control the time step.

BRIEF DESCRIPTION OF THE DRAWINGS

The structure and the technical means adopted by the present invention to achieve the above and other objects can be best understood by referring to the following detailed description of the preferred embodiments and the accompanying drawings, wherein:

FIG. 1a is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a neutral facial expression;

FIG. 1b is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a smile facial expression;

FIG. 1c is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a sad facial expression;

FIG. 1d is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a pout facial expression;

FIG. 1e is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a bitter smile facial expression;

FIG. 1f is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a pain facial expression;

FIG. 1g is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a wry facial expression;

FIG. 1h is one example of a facial expression for a person according to one embodiment of the present invention illustrating the sample case of a ferocious facial expression;

FIG. 2a is one example of a closed surface according to one embodiment of the present invention illustrating the sample case of a bimba closed surface;

FIG. 2b is one example of a closed surface according to one embodiment of the present invention illustrating the sample case of a cortex closed surface;

FIG. 2c is one example of a closed surface according to one embodiment of the present invention illustrating the sample case of a torso closed surface;

FIG. 2d is one example of a closed surface according to one embodiment of the present invention illustrating the sample case of a horse closed surface;

FIG. 3a shows convergence behavior of the computed harmonic energy produced by an explicit Euler method according to one exemplary embodiment of the present invention;

FIG. 3b shows convergence behavior of the computed harmonic energy produced by an implicit Euler method according to one exemplary embodiment of the present invention;

FIG. 3c shows convergence behavior of the computed harmonic energy produced by a quasi-implicit Euler method according to one exemplary embodiment of the present invention;

FIG. 4a shows convergence behavior of the computed harmonic energy produced by the phase-II quasi-implicit Euler method with an initial time step δt0 of 1 according to one exemplary embodiment of the present invention;

FIG. 4b shows convergence behavior of the computed harmonic energy produced by the phase-II quasi-implicit Euler method with an initial time step δt0 of 10 according to one exemplary embodiment of the present invention;

FIG. 4c shows convergence behavior of the computed harmonic energy produced by the phase-II quasi-implicit Euler method with an initial time step δt0 of 100 according to one exemplary embodiment of the present invention;

FIG. 4d shows convergence behavior of the computed harmonic energy produced by the phase-II quasi-implicit Euler method with an initial time step δt0 of 1000 according to one exemplary embodiment of the present invention;

FIG. 5a shows a human hand;

FIG. 5b shows the computed harmonic energies for various iterations of phase-II QIEM for the human hand according to one exemplary embodiment of the present invention;

FIG. 5c shows the computed time steps for various iterations of phase-II QIEM for the human hand according to one exemplary embodiment of the present invention;

FIG. 6a shows the computed harmonic energies for various iterations of the two-phase QIEM for the human torso in FIG. 2 c;

FIG. 6b shows the computed time steps for various iterations of the two-phase QIEM for the human torso in FIG. 2 c.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

In order to learn functions and features of the present invention, please refer to the following embodiments with detailed descriptions and the figures.

Conformal Mappings Spherical Conformal Mapping

First the spherical conformal mapping for genus zero closed surfaces from the point of view that a map is conformal if and only if it is harmonic is introduced. That means how to use the heat flow method to deform a mapping into the harmonic mapping under a special normalization condition is also introduced.

-   Suppose M is a triangular mesh of a genus zero closed surface with n     vertices {v1, . . . , vn}. All piecewise linear functions defined on     M is denoted by C^(PL)(M), which forms a linear space. -   Definition 1. (Discrete harmonic energy). Let f=(f₁, f₂, f₃): M→     ³ with f f₁, f₂, f₃ ∈ C^(PL)(M). The harmonic energy of f is defined     as

$\begin{matrix} {{ɛ_{h}(f)} = {\sum\limits_{l = 1}^{3}{ɛ_{h}\left( f_{l} \right)}}} & \left( {2a} \right) \\ {with} & \; \\ {{{ɛ_{h}\left( f_{l} \right)} = {\frac{1}{2}{\sum\limits_{{\lbrack{v_{i},v_{j}}\rbrack} \in M}{k_{ij}\left( {{f_{l}\left( v_{i} \right)} - {f_{l}\left( v_{j} \right)}} \right)}^{2}}}},{l = 1},2,3,} & \left( {2b} \right) \end{matrix}$

where {k_(ij)} forms a set of harmonic weights assigned on each edge [v_(i), v_(j)] ∈ M and is chosen such that the quadratic form of (2b) is positive definite.

-   Definition 2. Let f=(f₁, f₂, f₃): M→     ³ with f f₁, f₂, f₃ ∈ C^(PL)(M). The piecewise Laplacian operator     off is defined by

$\begin{matrix} {{\Delta_{d}f} = \left( {{\Delta_{d}f_{1}},{\Delta_{d}f_{2}},{\Delta_{d}f_{3}}} \right)} & (3) \\ {with} & \; \\ {{{\Delta_{d}{f_{l}\left( v_{i} \right)}} = {\sum\limits_{{\lbrack{v_{i},v_{j}}\rbrack} \in M}{k_{ij}\left( {{f_{l}\left( v_{i} \right)} - {f_{l}\left( v_{j} \right)}} \right)}}},{l = 1},2,3,} & \; \end{matrix}$

in which k_(ij) are the harmonic weights in (2b).

Let f (v) and n(f(v)) denote the image of the vertex v ∈ M and the normal on the target plane at f(v), respectively. Then the normal and tangent components of Δ_(d)f are defined as

(Δ_(d) f(υ))^(⊥)=<Δ_(d) f(υ), n(f(υ))>n(f(υ))   (4)

and

(Δ_(d) f(υ))^(∥)=Δ_(d) f(υ)−(Δ_(d) f(υ))^(⊥),   (5)

respectively, where <•, •> denotes the inner product in

³. Moreover, a map f: M₁→M₂ is harmonic, if and only if f only has a normal component, and the tangential component is zero, i.e.,

Δ_(d) f=(Δ_(d) f)^(⊥).   (6)

-   Remark 1. Note that if the following is taken

$\begin{matrix} {{a_{ii}:={\sum\limits_{{\lbrack{v_{i},v_{j}}\rbrack} \in M}k_{i,j}}},\mspace{14mu} {a_{ij}:={- k_{ij}}},{i \neq j},} & (7) \end{matrix}$

for i, j=1, n and define A≡[a_(ij)] ∈

^(n×n), then the discrete Laplacian operator Δ_(d)f in (2) can be represented as the matrix form

Δ_(d) f=Af=(Af ₁ ,Af ₂ ,Af ₃).   (8)

Many different ways are proposed to determine the edge weights k_(ij) in (2) so that the associated coefficient matrix A in (7) is symmetric positive semi-definite. A widely used edge weighting is the cotangent weighting. The matrix A associated with cotangent weights has been shown to be symmetric positive semi-definite.

A classical way to find the harmonic map f: M→

² is to minimize the discrete harmonic energy (2) by time evolution according to the nonlinear heat diffusion process

$\begin{matrix} {\frac{df}{dt} = {{- \Delta_{d}}{f.}}} & (9) \end{matrix}$

However, f(M, t) is constrained on the unit sphere

² so that it needs to project −Δ_(d)f onto the tangent plane of the sphere. Therefore, from (4)-(6), the evolution of f is according to the nonlinear heat diffusion equation:

$\begin{matrix} {\frac{{df}(v)}{dt} = {{- \left( {\Delta_{d}{f(v)}} \right)^{}} = {- {\left( {{\Delta_{d}{f(v)}} - {{\langle{{\Delta_{d}{f(v)}},{n\left( {f(v)} \right)}}\rangle}{n\left( {f(v)} \right)}}} \right).}}}} & (10) \end{matrix}$

One major difficulty is that the solution to the conformal mapping from M to

² is not unique but forms a Möbius group. In order to determine a unique solution, the additional zero mass-center constraint is required.

-   Definition 3. A mapping f: M₁→M₂ satisfies the zero mass-center     condition if and only if

∫_(M) ₁ fdσ_(M) ₁ =0,

where dσ_(M1) is the area element on M₁. All conformal maps from M to

² satisfying the zero mass-center constraint are unique up to the Euclidean rotation group.

Riemann Mapping

The spherical conformal mapping for genus zero closed surfaces can be utilized to find a conformal mapping (Riemann mapping) from a simply connected surface M with a single boundary ∂M to a two-dimensional (2D) unit disk

. The procedures for finding Riemann mapping are stated as follows. For a given simply connected triangular mesh M with boundary ∂M, there exists a symmetric closed surface Mc, called a double covering of M, which covers M twice, i.e., for each face in M, there are two preimages in Mc. Applying the spherical conformal mapping to Mc, a conformal map from Mc to the unit sphere

² is found. Next, a Möbius transformation τ:

→

${{\tau (z)} = \frac{{az} + b}{{cz} + d}},\mspace{11mu} a,b,c,{d \in {\mathbb{C}}},{{{ad} - {bc}} = 1},$

is used to adjust the conformal map such that ∂M is mapped to the equator of the unit sphere. Finally, the stereographic projection φ:

²→

${{\varphi \left( {x,y,z} \right)} = \left( {\frac{x}{1 - z},\frac{y}{1 - z}} \right)},{\left( {x,y,z} \right) \in ^{2}}$

is applied to map the lower hemisphere conformally to the unit disk

.

Quasi-Implicit Euler Method

Solving the steady state problems in (10) is the most time consuming step in finding conformal mappings. Thus an efficient algorithm is provided to solve (10).

The following definition is given for convenience.

-   Definition 1. Given

${{u \equiv \begin{bmatrix} u_{1} \\ \vdots \\ u_{n} \end{bmatrix}} \in {\mathbb{R}}^{n \times 3}},{{v \equiv \begin{bmatrix} v_{1} \\ \vdots \\ v_{n} \end{bmatrix}} \in {\mathbb{R}}^{n \times 3}},$

the operator

u, v

is denoted as

u, v

=diag(u ₁ v ₁ ^(T) , . . . , u _(n) v _(n) ^(T)).

Explicit and Implicit Euler Methods

For solving the steady state problem in (10), an explicit (forward) Euler method has been proposed by the following updating

$\begin{matrix} {\frac{f^{({m + 1})} - f^{(m)}}{\delta \; t} = {- \left( {{{{\Delta_{d}f^{(m)}} -} \prec {\Delta_{d}f^{(m)}}},{f^{(m)} \succ f^{(m)}}} \right)}} \\ {= {{- \left( {{{A -} \prec {A\; f^{(m)}}},{f^{(m)} \succ}} \right)}{f^{(m)}.}}} \end{matrix}$

Here the matrix A is the coefficient matrix of the Laplacian operator as in (7) with the edge weights k_(ij). The advantage of the explicit Euler method is that it only needs matrix-vector multiplications in each iteration. By choosing time step δt carefully, the associated energy can be monotonically diminished at each iteration, for example, when δt is chosen close to the square of the minimum of the edge lengths of M. However, in order to ensure numerical stability, the explicit technique always requires a very small time step which results in a significant drawback—a very slow convergence rate.

To remedy this obstacle, one may apply the implicit (backward) Euler method, which is A-stable over a wide range of time steps, to solve (10) as follows:

$\begin{matrix} \begin{matrix} {\frac{f^{({m + 1})} - f^{(m)}}{\delta \; t} = {- \left( {{{{\Delta_{d}f^{({m + 1})}} -} \prec {\Delta_{d}f^{({m + 1})}}},{f^{({m + 1})} \succ f^{({m + 1})}}} \right)}} \\ {{= {{- \left( {{{A -} \prec {A\; f^{({m + 1})}}},{f^{({m + 1})} \succ}} \right)}f^{({m + 1})}}},} \end{matrix} & (11) \end{matrix}$

or equivalently,

[I+(δt)A]f ^((m+1))−(δt)

Af ^((m+1)) , f ^((m+1))

f ^((m+1)) =f ^((m)).

The above equation may be rewritten as the nonlinear systems F(f^((m+1)), f^((m)))=0, where

F(f^((m+1)), f^((m)))≡−vec(f^((m)))+{(I₃

[I+(δt)A])−(δt)(I₃

Af^((m+1)), f^((m+1))

)}vec(f^((m+1))).   (12)

The unknown vectors f^((m+1)) of F(f^((m+1)), f^((m)) ) can be solved by Newton's method

vec(f _(i+1) ^((m+1)))=vec(f _(i) ^((m+1)))−J(f _(i) ^((m+1)))⁻¹ F(f _(i) ^((m+1)) , f ^((m)))   (13)

with f₀ ^((m+1))=f^((m)), where J (f^((m+1))) is the Jacobian matrix of F(f^((m+1)), f^((m))).

The implicit Euler method in (11) requires to solve a linear system if Newton's method is applied. Although the Jacobian matrix J (f_(i) ^((m+1))) can be reordered as a banded matrix so that direct methods can be applied, its size is enlarged to three times that of the matrix A.

Quasi-Implicit Euler Method

In order to avoid solving the nonlinear system in implicit Euler methods, some semi-implicit Euler methods have been proposed to improve the computational cost in each iteration. Based on these, the following quasi-implicit Euler method (QIEM) is proposed for solving the nonlinear heat diffusion equation in (10):

$\begin{matrix} {\frac{f^{({m + 1})} - f^{(m)}}{\delta \; t} = {- \left( {{{{\Delta_{d}f^{({m + 1})}} -} \prec {\Delta_{d}f^{(m)}}},{f^{(m)} \succ f^{({m + 1})}}} \right)}} \\ {= {{- \left( {{{A -} \prec {A\; f^{(m)}}},{f^{(m)} \succ}} \right)}{f^{({m + 1})}.}}} \end{matrix}$

That is, in each iteration, the new vector f^((m+1)) is generated by solving the linear system

[I+δt(A−

Af ^((m)) , f ^((m))

)]f ^((m+1)) =f ^((m)).   (14)

As an implicit Euler method, the QIEM has a wider range of stable time steps. Moreover, QIEM is much more efficient than the implicit Euler method by comparing the computational costs of solving the linear systems in (13) and (14).

Convergence Analysis of QIEM

In this section, the convergence of QIEM under the simplification of normalization of f^((m+1)) is analyzed.

The QIEM with simplification of the normalization of f^((m+1)) is stated as follows. Given f⁽⁰⁾ ∈

^(n×3) with

${{e_{i}^{T}f^{(0)}}}_{2} = \frac{1}{\sqrt{n}}$

for i=1, . . . , n, i.e., ∥f⁽⁰⁾∥_(F)=1, f^((m+1)) is defined by

$\begin{matrix} {{f^{({m + 1})} = {\frac{1}{\sqrt{n}}D_{m}^{{- 1}/2}A_{m}^{- 1}f^{(m)}}},} & (15) \end{matrix}$

for m=0, 1, . . . , where

A _(m) =I+δt(A−

Af ^((m)) , f ^((m))

),   (16)

D _(m) =

A _(m) ⁻¹ f ^((m)) , A _(m) ⁻¹ f ^((m))

.   (17)

-   Note that

${{e_{i}^{T}f^{({m + 1})}}}_{2} = \frac{1}{\sqrt{n}}$

for i=1, . . . , n and ∥f ^((m+1))∥_(F)=1.

Let us consider the Schur decomposition

Q _(m)Λ_(m) Q _(m) ^(T) =A−

Af ^((m)) , f ^((m))

,   (18)

where Q_(m)≡[_(q1,m, . . . , qn,m)]^(T) is orthonormal and Λ_(m)=diag(λ₁ ^((m)), . . . , λ_(n) ^((m))) with λ_(i) ^((m)) being the eigenvalues. It is assumed λ_(i) ^((m))≠0 for i=1, . . . , n. Then

$\begin{matrix} \begin{matrix} {{\frac{1}{\sqrt{n}}{D_{m}^{{- 1}/2}}_{2}} = {\frac{1}{\sqrt{n}}{\max\limits_{i}\left\{ {{e_{i}^{T}A_{m}^{- 1}f^{(m)}}}_{2}^{- 1} \right\}}}} \\ {= {\frac{1}{\sqrt{n}}{\max\limits_{i}\left\{ {{{q_{i,m}^{T}\left( {I + {\delta \; t\; \Lambda_{m}}} \right)}^{- 1}Q_{m}^{T}f^{(m)}}}_{2}^{- 1} \right\}}}} \\ {{= {\frac{\delta \; t}{\sqrt{n}}{\max\limits_{i}\left\{ {{{q_{i,m}^{T}\; \Lambda_{m}^{- 1}Q_{m}^{T}f^{(m)}} + {O_{m}\left( \frac{1}{\delta \; t} \right)}}}_{2}^{- 1} \right\}}}},} \end{matrix} & (19) \end{matrix}$

where Om (1/δt) is dependent on (i,m). Given a small positive value η_(m)>0, there exists T₀>0 such that

$\begin{matrix} {{{{O_{m}\left( \frac{1}{\delta \; t} \right)}}_{2} \leq \eta_{m}},{\forall{{\delta \; t} \geq T_{0}}}} & (20) \end{matrix}$

which implies that

$\begin{matrix} {{\frac{1}{\sqrt{n}}{D_{m}^{{- 1}/2}}_{2}} \leq {\frac{\delta \; t}{\sqrt{n}}{\max\limits_{i}\left\{ {{{q_{i}^{T}\; \Lambda_{m}^{- 1}Q_{m}^{T}f^{(m)}}}_{2 - \eta_{m}}}^{- 1} \right\}}} \equiv {\frac{\delta \; t}{\sqrt{n}}{\frac{1}{a_{m}}.}}} & (21) \end{matrix}$

By the definition of f^((m+1)) in (15), it holds that

$\begin{matrix} {\begin{matrix} {E_{m + 1} \equiv {f^{({m + 1})} - f^{(m)}}} \\ {= {{\frac{1}{\sqrt{n}}D_{m}^{{- 1}/2}A_{m}^{- 1}E_{m}} + {\frac{1}{\sqrt{n}}{D_{m}^{{- 1}/2}\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)}f^{({m - 1})}} +}} \\ {{\frac{1}{\sqrt{n}}\left( {D_{m}^{{- 1}/2} - D_{m - 1}^{{- 1}/2}} \right)A_{m - 1}^{- 1}{f^{({m - 1})}.}}} \end{matrix}{Assume}} & (22) \\ {{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} = {{\min\limits_{i}{{1 + {\delta \; t\; \lambda_{i}^{(m)}}}}} > 1}} & (23) \end{matrix}$

for δt>0, where λ _(m)>0.

-   Lemma 1. Let A_(m), D_(m) and E_(m) be defined in (16), (17) and     (22), respectively. Let a_(m) and λ _(m) be defined in (21) and     (23), respectively. Assume

$\begin{matrix} {{{{\sqrt{n}a_{m}{\underset{\_}{\lambda}}_{m}} - 3} > 0.}{Then}{{\frac{1}{\sqrt{n}}{{D_{m}^{{- 1}/2}A_{m}^{- 1}E_{m}}}_{F}} < {\frac{1}{3}{E_{m}}_{F}}}} & (24) \end{matrix}$

for all

${\delta \; t} > T_{1} \equiv {\max \left\{ {\frac{\sqrt{n}a_{m}}{{\sqrt{n}a_{m}{\underset{\_}{\lambda}}_{m}} - 3},T_{0}} \right\}}$

where T0 is defined in (20).

-   Proof. From the assumption in (24), the following is obtained

(√{square root over (n)}a_(m) λ _(m)−3)δt>√{square root over (n)}a _(m), for δt>T ₁

which is equivalent to

−3δt>√{square root over (n)}a _(m)(1−δtλ _(m))

and then

3δt<√{square root over (n)}a _(m)|1−δtλ _(m)|.   (25)

Combing the results in (21) and (25), it holds that

${\frac{1}{\sqrt{n}}{{D_{m}^{{- 1}/2}A_{m}^{- 1}E_{m}}}_{F}} \leq {\frac{1}{\sqrt{n}}\frac{\delta \; t}{a_{m}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}}}{E_{m}}_{F}} < {\frac{1}{3}{{E_{m}}_{F}.}}$

-   Lemma 2. Assume that

α_(m) ≡a _(m) √{square root over (n)}λ _(m) λ _(m−1)−3b(1+1/√{square root over (n)})>0,   (26)

where b≡max_(1≦i≦n) ∥e _(i) ^(T) A∥_(2·). Then

$\begin{matrix} {{\frac{1}{\sqrt{n}}{{{D_{m}^{{- 1}/2}\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)}f^{({m - 1})}}}_{F}} < {\frac{1}{3}{E_{m}}_{F}}} & (27) \end{matrix}$

for all

${{\delta \; t} > T_{2} \equiv {\max \left\{ {T_{0},\frac{{- \beta_{m}} + \sqrt{\beta_{m}^{2} - {4\; \alpha_{m}\gamma_{m}}}}{2\; \alpha_{m}}} \right\}}},$

where T₀ is defined in (20) and

β_(m) =−a _(m) √{square root over (n)}(λ _(m)+λ _(m−1)), γ_(m)=α_(m) √{square root over (n)}.

Proof. By the definition of A_(m) in (16), it holds that

$\begin{matrix} {\mspace{45mu} {\begin{matrix} {\mspace{11mu} {\begin{matrix} {\frac{1}{\sqrt{n}}D_{m}^{{- 1}/2}} \\ {\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)f^{({m - 1})}} \end{matrix} = {\frac{1}{\sqrt{n}}D_{m}^{{- 1}/2}{A_{m}^{- 1}\left( {A_{m - 1} - A_{m}} \right)}}}} \\ {{A_{m - 1}^{- 1}f^{({m - 1})}}} \\ {= {\frac{\delta \; t}{\sqrt{n}}D_{m}^{{- 1}/2}{A_{m}^{- 1}\left( {{\prec {A\; f^{(m)}}},} \right.}}} \\ \left. {{f^{(m)} \succ {- {\prec {A\; f^{({m - 1})}}}}},{f^{({m - 1})} \succ}} \right) \\ {{A_{m - 1}^{- 1}f^{({m - 1})}}} \\ {= {\frac{\delta \; t}{\sqrt{n}}D_{m}^{{- 1}/2}{A_{m}^{- 1}\left( {{\prec {A\; f^{(m)}}},} \right.}}} \\ \left. {{E_{m} \succ {+ {\prec {AE}_{m}}}},{f^{({m - 1})} \succ}} \right) \\ {{A_{m - 1}^{- 1}{f^{({m - 1})}.}}} \end{matrix}\mspace{79mu} {Since}\begin{matrix} {\mspace{79mu} {\begin{matrix} {{{{\prec {A\; f^{(m)}}},{E_{m} \succ}}}_{2} \leq} \\ {\max\limits_{i}{{{e_{i}^{T}A}}_{2}{f^{(m)}}_{F}{E_{m}}_{F}}} \end{matrix} = {\max\limits_{i}{{e_{i}^{T}A}}_{2}}}} \\ {{E_{m}}_{F}} \\ {\equiv {b{E_{m}}_{F}}} \end{matrix}\mspace{20mu} {and}\begin{matrix} {\mspace{79mu} {{{{\prec {AE}_{m}},{f^{({m - 1})} \succ}}}_{2} = {{\max\limits_{i}{{\left( {e_{i}^{T}{AE}_{m}} \right)\left( {e_{i}^{T}f^{({m - 1})}} \right)^{\prime}}}} \leq}}} \\ {{{\max\limits_{i}{{{e_{i}^{T}{AE}_{m}}}_{2}{\max\limits_{i}{{e_{i}^{T}f^{({m - 1})}}}_{2}}}} \leq}} \\ {{\max\limits_{i}{{{e_{i}^{T}A}}_{2}\frac{1}{\sqrt{n}}{E_{m}}_{F}}}} \\ {{\equiv {\frac{b}{\sqrt{n}}{E_{m}}_{F}}},} \end{matrix}\mspace{20mu} {thus}{{{\frac{1}{\sqrt{n}}{{{D_{m}^{{- 1}/2}\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)}f^{({m - 1})}}}_{F}} \leq {\frac{\left( {\delta \; t} \right)^{2}}{\sqrt{n}}\frac{b + {b/\sqrt{n}}}{{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}a_{m}}{E_{m}}_{F}{f^{({m - 1})}}_{F}}} = {\frac{\left( {\delta \; t} \right)^{2}}{\sqrt{n}}\frac{b + {b/\sqrt{n}}}{{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}a_{m}}{{E_{m}}_{F}.}}}}} & (28) \end{matrix}$

-   By the definitions of α_(m), β_(m) and γ_(m), it follows that

β_(m) ²4α_(m)γ_(m) =na _(m) ²(λ _(m)−λ _(m−1))²+12√{square root over (n)}a _(m) b(1+1/√{square root over (n)})>0.

-   Using the assumption αm>0, for all δt>T2, the following is obtained

${\left( {{\delta \; t} - \frac{{- \beta_{m}} + \sqrt{\beta_{m}^{2} - {4\; \alpha_{m}\gamma_{m}}}}{2\; \alpha_{m}}} \right)\mspace{14mu} \left( {{\delta \; t} - \frac{{- \beta_{m}} - \sqrt{\beta_{m}^{2} - {4\; \alpha_{m}\gamma_{m}}}}{2\; \alpha_{m}}} \right)} > 0$

which implies that

√{square root over (n)}a _(m) [1−(λ _(m)+λ _(m−1))δt+λ _(m) λ _(m−1)(δt)²]>3b(1+1/√{square root over (n)})(δt)².

-   Consequently,

$\begin{matrix} {{\frac{\left( {\delta \; t} \right)^{2}}{\sqrt{n}}\frac{b + {b/\sqrt{n}}}{{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}a_{m}}} < {\frac{1}{3}.}} & (29) \end{matrix}$

-   Substituting (29) into (28), the result in (27) is obtained. -   Lemma 3. Assume that

$\begin{matrix} {{{\min\limits_{1 \leq i \leq n}{{q_{i,m}^{T}A_{m}^{- 1}Q_{m}^{T}f^{(m)}}}_{2}} > \eta_{m}},} & (30) \end{matrix}$ 3(|λ _(m)|+|λ _(m−1)|)(b+b/√ n )+6|λ _(m) λ _(m−1) |<√{square root over (n)}a _(m) a _(m−1) c _(m) λ _(m) ² λ _(m−1) ²|λ _(m)|,   (31)

where η_(m) is defined in (20) and

$c_{m} = {\min\limits_{i}{\left( {{{q_{i,m}^{T}A_{m}^{- 1}Q_{m}^{T}f^{(m)}}}_{2} + {{q_{i,{m - 1}}^{T}A_{m - 1}^{- 1}Q_{m - 1}^{T}f^{({m - 1})}}}_{2} - \eta_{m} - \eta_{m - 1}} \right).}}$

-   Then there exists T₃>0 such that

${\frac{1}{\sqrt{\eta}}{{\left( {D_{m}^{{- 1}/2} - D_{m - 1}^{{- 1}/2}} \right)A_{m - 1}^{- 1}f^{({m - 1})}}}_{F}} < {\frac{1}{3}{E_{m}}_{F}}$

for δt≧T₃.

-   Proof. The third term on the right hand side of (22) is rewritten as

$\begin{matrix} {{\frac{1}{\sqrt{\eta}}\left( {D_{m}^{{- 1}/2} - D_{m - 1}^{{- 1}/2}} \right)A_{m - 1}^{- 1}f^{({m - 1})}} = {\frac{1}{\sqrt{\eta}}\left( {D_{m} - D_{m - 1}} \right)\left\{ {D_{m}^{{- 1}/2}{D_{m - 1}^{{- 1}/2}\left( {D_{m}^{1/2} + D_{m - 1}^{1/2}} \right)}^{- 1}} \right\} A_{m - 1}^{- 1}{f^{({m - 1})}.}}} & (32) \end{matrix}$

-   From (19) and (20), the following is obtained

w i ≡  e i T  ( D m 1 / 2 + D m - 1 1 / 2 )  e i =  1 δ   t  {  q i , m T  Λ m - 1  Q m T  f ( m ) + m  ( 1 δ   t )  2 +   q i , m - 1 T  Λ m - 1 - 1  Q m - 1 T  f ( m - 1 ) + m - 1  ( 1 δ   t )  2 } ≥  1 δ   t  {  q i , m T  Λ m - 1  Q m T  f ( m )  2 +  q i , m - 1 T  Λ m - 1 - 1  Q m - 1 T  f ( m - 1 )  2 -  η m - η m - 1 }

for δt≧T₀, which implies that

$\begin{matrix} \begin{matrix} {{\left( {D_{m}^{1/2} + D_{m - 1}^{1/2}} \right)^{- 1}}_{2} = {\max\limits_{1 \leq i \leq n}w_{i}^{- 1}}} \\ {= {\left\{ {\max\limits_{1 \leq i \leq n}w_{i}} \right\}^{- 1} \leq {\delta \; t\left\{ \min\limits_{i} \right.}}} \\ {\left( {{{q_{i,m}^{T}\Lambda_{m}^{- 1}Q_{m}^{T}f^{(m)}}}_{2} +} \right.} \\ {{{{q_{i,{m - 1}}^{T}\Lambda_{m - 1}^{- 1}Q_{m - 1}^{T}f^{({m - 1})}}}_{2} - \eta_{m} -}} \\ \left. \left. \eta_{m - 1} \right) \right\}^{- 1} \\ {\equiv {\frac{\delta \; t}{c_{m}}.}} \end{matrix} & (33) \end{matrix}$

-   On the other hand,

$\begin{matrix} \begin{matrix} {{{D_{m} - D_{m - 1}} = {\prec {A_{m}^{- 1}f^{(m)}}}},{{A_{m}^{- 1}f^{(m)}} \succ {- {\prec {A_{m}^{- 1}f^{(m)}}}}},{{A_{m}^{- 1}f^{(m)}} \succ}} \\ {{= {\prec {A_{m}^{- 1}f^{(m)}}}},{{\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)f^{(m)}} \succ +}} \\ {{{\prec {A_{m}^{- 1}f^{(m)}}},{{A_{m - 1}^{- 1}E_{m}} \succ +}}} \\ {{{\prec {A_{m}^{- 1}E_{m}}},{{A_{m - 1}^{- 1}f^{({m - 1})}} \succ +}}} \\ {{{\prec {\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)f^{({m - 1})}}},{{A_{m - 1}^{- 1}f^{({m - 1})}} \succ .}}} \end{matrix} & (34) \end{matrix}$

-   From (28), it follows that

$\begin{matrix} {{{A_{m}^{- 1} - A_{m - 1}^{- 1}}}_{2} \leq {\frac{\left( {b + {b/\sqrt{n}}} \right)\delta \; t}{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}{{E_{m}}_{F}.}}} & (35) \end{matrix}$

-   From (23), (34) and (35) with the results, the following is obtained

$\begin{matrix} {{{D_{m} - D_{m - 1}}}_{2} \leq {{{{\prec {A_{m}^{- 1}f^{(m)}}},{{\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)f^{(m)}} \succ}}}_{2} + {{{\prec {A_{m}^{- 1}f^{(m)}}},{{A_{m - 1}^{- 1}E_{m}} \succ}}}_{2} + {{{\prec {A_{m}^{- 1}E_{m}}},{{A_{m - 1}^{- 1}f^{({m - 1})}} \succ}}}_{2} + {{{\prec {\left( {A_{m}^{- 1} - A_{m - 1}^{- 1}} \right)f^{({m - 1})}}},{{A_{m - 1}^{- 1}f^{({m - 1})}} \succ}}}_{2}} \leq {{\left( {{A_{m}^{- 1}}_{2} + {A_{m - 1}^{- 1}}_{2}} \right){{A_{m}^{- 1} - A_{m - 1}^{- 1}}}_{2}} + {2{A_{m}^{- 1}}_{2}{A_{m - 1}^{- 1}}_{2}{E_{m}}_{2}}} \leq {{\left( {\frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} + \frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}} \right)\frac{\left( {b + {b/\sqrt{n}}} \right)\delta \; t}{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}{E_{m}}_{F}} + {\frac{2}{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}{{E_{m}}_{F}.}}}} & (36) \end{matrix}$

-   Using (23), (34) and (35) and the assumption (31), (32) implies that

${{{\frac{1}{\sqrt{n}}{{\left( {D_{m}^{{- 1}/2} - D_{m - 1}^{{- 1}/2}} \right)A_{m - 1}^{- 1}f^{({m - 1})}}}_{F}} \leq {{{D_{m} - D_{m - 1}}}_{2}{D_{m}^{{- 1}/2}}_{2}{D_{m - 1}^{{- 1}/2}}_{2}{\left( {D_{m}^{1/2} + D_{m - 1}^{1/2}} \right)^{- 1}}_{2}{A_{m - 1}^{- 1}}_{2}} \leq {\frac{1}{\sqrt{n}}\frac{\delta \; t}{a_{m}}\frac{\delta \; t}{a_{m - 1}}\frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}}\frac{\delta \; t}{c_{m}}\left\{ {\frac{2}{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}} + {\left( {\frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} + \frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}} \right)\frac{\left( {b + {b/\sqrt{n}}} \right)\delta \; t}{{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} \cdot {{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}}}} \right\} {E_{m}}_{F}}} = {{\frac{\left( {\delta \; t} \right)^{2}}{\sqrt{n}a_{m}a_{m - 1}{c_{m}\left( {1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}} \right)}^{2}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}} \times \left\{ {2 + {\left( {\frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m}}}} + \frac{1}{{1 - {\delta \; t{\underset{\_}{\lambda}}_{m - 1}}}}} \right)\left( {b + {b/\sqrt{n}}} \right)\delta \; t}} \right\}} < {\frac{1}{3}{E_{m}}_{F}}}},\mspace{20mu} \left( {{by}\mspace{14mu} {the}\mspace{14mu} {assumption}\mspace{14mu} (31)} \right)$

for δt≧T₃ with large enough T₃.

-   Using the results in Lemmas 1, 2 and 3 into (22), the decrease of     the error ∥E_(m)∥_(F) is easily shown. -   Theorem 1. Assume that inequalities of (24), (26), (30) and (31)     hold. -   Then there exists T>0 such that

∥E _(m)∥_(F) ≡∥f ^((m+1)) −f ^((m))∥_(F) <∥f ^((m)) −f ^((m−1))∥_(F) ≡∥E _(m−1)∥_(F)

for all δt>T.

-   Now, a more simplified case for f^((m)) and f^((m+1)) without     normalization is considered. Assume

∥e _(i) ^(T) f ^((m−1))∥₂=1, for i=1, . . . , n.   (37)

Let

f ^((m)) =A _(m−1) ⁻¹ f ^((m−1)) , f ^((m+1)) =A _(m) ⁻¹ f ^((m)),

where A_(m) is defined in (16) and

α_(m)=λ _(m−1) λ _(m),   (38)

β_(m)=λ _(m)+λ _(m−1)+4b√{square root over (n)}.   (39)

where λ _(m) is defined in (23). Then f^((m−1)), f^((m)) and f^((m+1)) satisfy the following result.

-   Theorem 2. Assume α_(m), β_(m) ²−4α_(m)>0. If δt satisfies that

$\begin{matrix} {{\delta \; t} > \left\{ {\begin{matrix} {{\left( {\sqrt{n} + 1} \right)/{\underset{\_}{\lambda}}_{m}},} & {{{{if}\mspace{14mu} {\underset{\_}{\lambda}}_{m}} > 0};} \\ {{\left( {\sqrt{n} - 1} \right)/\left( {- {\underset{\_}{\lambda}}_{m}} \right)},} & {{{{if}\mspace{14mu} {\underset{\_}{\lambda}}_{m}} < 0},} \end{matrix}{and}} \right.} & (40) \\ {{{{\delta \; t} > \frac{\beta_{m} + \sqrt{\beta_{m}^{2} - {4\; \alpha_{m}}}}{4\; \alpha_{m}}},{then}}{{{f^{({m + 1})} - f^{(m)}}}_{F} < {{{f^{(m)} - f^{({m - 1})}}}_{F}.}}} & (41) \end{matrix}$

-   Proof. From (16), the following is obtained

$\begin{matrix} \begin{matrix} {E_{m + 1} \equiv {f^{({m + 1})} - f^{(m)}}} \\ {= {{A_{m}^{- 1}f^{(m)}} - {A_{m - 1}^{- 1}f^{({m - 1})}}}} \\ {= {{A_{m}^{- 1}f^{(m)}} - {A_{m}^{- 1}f^{({m - 1})}} + {A_{m}^{- 1}f^{({m - 1})}} - {A_{m - 1}^{- 1}f^{({m - 1})}}}} \\ {= {{A_{m}^{- 1}E_{m}} + {{A_{m}^{- 1}\left( {A_{m - 1} - A_{m}} \right)}A_{m - 1}^{- 1}f^{({m - 1})}}}} \\ {= {{A_{m}^{- 1}E_{m}} + {\left( {\delta \; t} \right){A_{m}^{- 1}\left( {{\prec {A\; f^{(m)}}},{f^{(m)} \succ {- {\prec {A\; f^{({m - 1})}}}}},} \right.}}}} \\ {\left. {f^{({m - 1})} \succ} \right)A_{m - 1}^{- 1}f^{({m - 1})}} \\ {= {{A_{m}^{- 1}E_{m}} + {\left( {\delta \; t} \right){A_{m}^{- 1}\left( {{\prec {A\; f^{(m)}}},{f^{(m)} \succ {- {\prec {A\; f^{(m)}}}}},} \right.}}}} \\ {{{f^{({m - 1})} \succ {+ {\prec {A\; f^{(m)}}}}},{f^{({m - 1})} \succ {- {\prec {A\; f^{({m - 1})}}}}},}} \\ {\left. {f^{({m - 1})} \succ} \right)A_{m - 1}^{- 1}{f^{({m - 1})}.}} \end{matrix} & (42) \end{matrix}$

-   By QIEM definition 1, it holds that

$\begin{matrix} {{{{\prec {A\; f^{(m)}}},{f^{({m - 1})} \succ {- {\prec {A\; f^{({m - 1})}}}}},{f^{({m - 1})} \succ}}}_{2} = {{{\prec {AE}_{m}},{f^{({m - 1})} \succ}}}_{2}} \\ {= {\max\limits_{1 \leq i \leq n}{{e_{i}^{T}{AE}_{m}}}}} \\ {{{\left( f^{({m - 1})} \right)^{T}e_{i}}}.} \end{matrix}$

-   By the assumption in (37), the following is obtained

$\begin{matrix} {{{{\prec {A\; f^{(m)}}},{f^{({m - 1})} \succ {- {\prec {A\; f^{({m - 1})}}}}},{f^{({m - 1})} \succ}}}_{2} \leq {\max\limits_{1 \leq i \leq n}{{e_{i}^{T}{AE}_{m}}}_{2}} \leq {\max\limits_{1 \leq i \leq n}{{{e_{i}^{T}A}}_{2}{E_{m}}_{F}}} \equiv {b{{E_{m}}_{F}.}}} & (43) \end{matrix}$

-   On the other hand, using f(m) and (16), (23) and (37), it holds that

$\begin{matrix} {{{{\prec {A\; f^{(m)}}},{f^{(m)} \succ {- {\prec {A\; f^{(m)}}}}},{f^{({m - 1})} \succ}}}_{2} = {{\max\limits_{1 \leq i \leq n}{{e_{i}^{T}{AA}_{m}^{- 1}f^{({m - 1})}E_{m}^{T}e_{i}}}} \leq {\max\limits_{1 \leq i \leq n}{{{e_{i}^{T}A}}_{2}{A_{m}^{- 1}}_{2}{f^{({m - 1})}}_{2}{E_{m}}_{F}}} \leq {\frac{b\sqrt{n}}{{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}}{{E_{m}}_{F}.}}}} & (44) \end{matrix}$

-   The assumption in (40) implies that

$\begin{matrix} {{\frac{\sqrt{n}}{{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}} < {1\mspace{14mu} {and}\mspace{14mu} \frac{1}{{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}}} < \frac{1}{\sqrt{n}} < \frac{1}{2}},{{{if}\mspace{14mu} n} > 4.}} & (45) \end{matrix}$

-   Consequently, from (44),

∥

Af ^((m)) , f ^((m))

−

Af ^((m)) , f ^((m−1))

∥₂ <b∥E _(m)∥_(F).   (46)

-   From (42), (43) and (46), the following is obtained

$\begin{matrix} {{E_{m + 1}}_{F} < {{{A_{m}^{- 1}}_{2}{E_{m}}_{F}} + {2\; b\mspace{11mu} \left( {\delta \; t} \right){A_{m}^{- 1}}_{2}{A_{m - 1}^{- 1}}_{2}{f^{({m - 1})}}_{2}{E_{m}}_{F}}} < {\left\{ {\frac{1}{{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}} + \frac{2\; b\mspace{11mu} \left( {\delta \; t} \right)\mspace{11mu} \sqrt{n}}{{{1 - {{\underset{\_}{\lambda}}_{m - 1}\delta \; t}}}\mspace{14mu} {{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}}}} \right\} {{E_{m}}_{F}.}}} & (47) \end{matrix}$

-   By the assumption that α_(m), β_(m) ²−4α_(m)>0, for all

${{\delta \; t} > \frac{\beta_{m} + \sqrt{\beta_{m}^{2} - {4\alpha_{m}}}}{2\alpha_{m}}},{{\left( {{\delta \; t} - \frac{\beta_{m} + \sqrt{\beta_{m}^{2} - {4\alpha_{m}}}}{2\alpha_{m}}} \right)\mspace{14mu} \left( {{\delta \; t} - \frac{\beta_{m} - \sqrt{\beta_{m}^{2} - {4\alpha_{m}}}}{2\alpha_{m}}} \right)} > 0}$

which implies that

$\begin{matrix} {{\left\lbrack {1 - {\left( {{\underset{\_}{\lambda}}_{m} + {\underset{\_}{\lambda}}_{m - 1}} \right)\mspace{11mu} \delta \; t} + {{\underset{\_}{\lambda}}_{m}{\underset{\_}{\lambda}}_{m - 1}\; \left( {\delta \; t} \right)^{2}}} \right\rbrack > {4\; b\sqrt{n}\delta \; t}}{And}{\frac{2\; b\mspace{11mu} (\; {\delta t})\mspace{11mu} \sqrt{n}}{{{1 - {{\underset{\_}{\lambda}}_{m - 1}\delta \; t}}}\mspace{14mu} {{1 - {{\underset{\_}{\lambda}}_{m}\delta \; t}}}} < {\frac{1}{2}.}}} & (48) \end{matrix}$

-   Substituting (45) and (48) into (47), it holds that     ∥E_(m+1)∥_(F)<∥E_(m)∥_(F).

Practical Implementation

How to efficiently compute the steady state solution of (10) with zero mass-center normalization by QIEM will be discussed.

1. The Initial Mapping f⁽⁰⁾

For the spherical conformal mapping on genus-zero closed surfaces, the initial map f⁽⁰⁾ in (14) is to construct an one-to-one and onto smooth map from a given genus-zero closed surface to the unit sphere. In practice, the Gauss map is used as the initial map f⁽⁰⁾, which is defined as follows:

-   Definition of Gauss map.     : M→     ², G(v)=n(v), where n(v) is the unit normal vector at v ∈ M.     The corresponding Gauss map can be computed, for instance, as the     weighted sum of normals on the adjacent faces weighed by their     areas.

For the case of simply connected surfaces with a single boundary, i.e., computing the Riemann mapping, the idea described in Connectivity shapes in proceedings of IEEE Visualization 2001, pp. 135-142, IEEE Computer Society, 2001 to construct a harmonic-type initial map. Specifically, first the open surface is parametrized onto a unit disk by the harmonic map. Using the stereographic projection, the resulting mesh is mapped to a lower hemisphere. Then the harmonic-type initial map is obtained by reflecting the hemisphere's image along the equatorial plane to built a full unit sphere.

2. Adaptive Time Step Control

The convergence of solving the steady state solution of the time-dependent differential equations with fixed time step δt may be very slow. In order to accelerate the convergence, a time step controlling scheme is proposed to control δt in (14) for each iteration m. As a consequence, the sequence {f^((m))} can be convergent as soon as possible.

Let δt₀ be an initial time step and δt_(max) denote the maximal time step so that δt is set to be δt_(max) if δt>δ_(max). The strategy of choosing δt in each iteration is based on the decrement of the harmonic energy ε_(h)(f), i.e., the time steps are chosen such that the harmonic energy is always decreasing in each iteration. This objective can be achieved by the descending and accelerating strategies of time step controls. The descending strategy is used to determine δt so that ε_(h)(f^((m+1))) is always decreasing, i.e., ε_(h)(f^((m+1)))−ε_(h)(f^((m)))<ε, in each m=0, 1, . . . , k for some k and tolerance ε. Given ε, δt₀ and an increment α, the strategy is described as follows.

-   (S₂.a) If ε_(h)(f^((m+1)))−ε_(h)(f^((m)))≧ε, then introduce δt to be     δt:=1/2δt and recompute f^((m+1)) from (14). This step prevents the     sequence {f^((m))} converges to a map whose harmonic energy is a     local extreme. -   (S₂.b) If the new ε_(h)(f^((m+1))) is still larger than     ε_(h)(f^((m))) when δt has been repeatedly reduced three times, then     the initial δt₀ is introduced to be δt₀:=δt₀/α and the algorithm is     restarted with m=0.

The iteration in (14) with the descending strategy is repeated till ε_(h)(f^((m+1))) is closed enough to ε_(h)(f^((m))). After the descending strategy, the following accelerating strategy is used to increase δt and to speed up the convergence further.

-   (S₂.c) If Εh(f(m+1)) satisfies that     |ε_(h)(f^((m+1)))−ε_(h)(f^((m)))|<10×|ε_(h)(f^((m)))−ε_(h)(f^((m−1)))|,     then δt is increased by δt:=min(δt_(max), δt×α); otherwise, δt is     too large, will be reduced by δt:=δt/α and f^((m+1)) is recomputed     from (14).

3. Two-Phase Quasi-Implicit Euler Method

To determine the initial map f⁽⁰⁾ in (14) is a crucial cornerstone to implement the QIEM for the evolution of conformal mappings on the unit sphere or on the unit disk. If f⁽⁰⁾ is not close to the steady state solution, then it is difficult to find a large stable convergence region of the time steps for solving (10) by QIEM with time step controlling strategies (S₂.a), (S₂.b) and (S₂.c). To remedy this drawback, a two-phase QIEM is proposed. If f⁽⁰⁾ is not close to the steady state solution, the heuristic phase-I QIEM

${\frac{{f^{({m + 1})} - f^{(m)}}}{\delta \; t} = {{{- \Delta_{d}}f^{({m + 1})}} = {- {Af}^{\; {({m + 1})}}}}},{i.e.},{{\left( {I + {\left( {\delta \; t} \right)A}} \right)f^{({m + 1})}} = f^{(m)}}$

is used to compute f^((m+1)). When f^((m)) is close to the steady state solution, switch to (14) with strategies (S₂.a), (S₂.b) and (S₂.c) for computing f^((m+1)). This is called the phase-II QIEM.

To set up a robust phase-I QIEM, an adaptive method is also proposed, just as the phase-II QIEM, to control δt in each iteration. The strategy described as follows is repeated until the difference |ε_(h)(f^((m+1)))−ε_(h)(f^((m)))| of the energy is less than a given tolerance ε₃.

-   (S₁.a) By the definition of A in (7), the row sums of A are equal to     zero which implies that there is a trivial solution f such that (A−     Af, f     )f=0 and ε_(h)(f)=0. To avoid f^((m)) convergent to this trivial     solution, reset δt₀ to be the current δt and restart the algorithm     with the original f⁽⁰⁾ and new δt₀ when ε_(h)(f^((m))) is less than     a given small tolerance ε₂. -   (S₁.b) If ε_(h)(f^((m+1))) does not decrease, then δt is reduced to     be δt:=max(1, 1/2δt) and recompute f^((m+1)) from (14). -   (S₁.c) If the difference of the energy is still larger than or equal     to ε₃ when m has been larger than a given maximal iteration number     m_(max), reset δt₀ to be the current δt and restart the algorithm     with original initial closed surface f⁽⁰⁾ and new δt₀.

It is important to determine δt₀ for solving the steady state problem. In general, the explicit Euler method has extremely narrow stability region of time steps which leads to the extremely slow convergence. The implicit Euler method has considerably wider stability region. However, it involves nonlinear systems that have to be solved by iterative solvers. If the time step is larger, then the initial guess of the iterative solver must be unrealistically close to the solution of the nonlinear systems. For two phases of QIEM, a heuristic method is proposed in Algorithm 1 to determine an initial time step. Comparing the initial time step in Algorithm 1 with those for the explicit and implicit Euler methods, the proposed time step is large.

-   The algorithm of “Two-phase quasi-implicit Euler method” is stated     as the following Algorithm 1. -   Input: A genus zero triangular mesh M, an initial map f⁽⁰⁾ with the     harmonic energy ε_(h)(f⁽⁰⁾) and a threshold δΕ for the energy     difference. -   Output: Convergent steady state solution f^((m)) with zero mass     center.     -   1: Compute the areas         _(i), i=1, . . . , n, of all faces in M;     -   2: Resort         _(i) for i=1, . . . , n;     -   3: Compute the average         of the areas         _(i) for i=0.45 n, . . . , 0.55 n;     -   4: Find the minimal positive integer k such that         ^(k/3-7)·≦3500.     -   5: Compute δt₀=         ^(k/2-7)·.     -   6: if f⁽⁰⁾ is not close to the steady state solution then     -   7: repeat     -   8: Solve the linear system

(I+(δt)A)f ^((m+1)) =f ^((m)).

-   -   9: Compute the mass sphere center c:

$\begin{matrix} {{c = {{\left( {\sum\limits_{i = 1}^{n}\; {\left( v_{i} \right)}} \right)^{- 1}\left\lbrack {\left( v_{i} \right)\mspace{14mu} \ldots \mspace{14mu} \left( v_{n} \right)} \right\rbrack}f^{({m + 1})}}},} & (49) \end{matrix}$

-   -   -   where             (υ_(i)) is the sum of areas for the faces which have the             common vertex v_(i), i=1, . . . , n;

    -   10: Normalize f^((m+1)) by (f^((m+1))−c)/∥f^((m+1))−c∥ such that         the mass center of the resulting f^((m+1)) is at the unit sphere         center;

    -   11: Compute δt according to strategies (S₁.a), (S₁.b) and (S₁.c)         with m_(max)=10, ε₂=0.5 and ε₃=0.1;

    -   12: until |ε_(h)(f^((m+1)))−ε_(h)(f^((m)))|<0.1

    -   13: Set f⁽⁰⁾=f^((m+1)) and m=0;

    -   14: end if

    -   15: repeat

    -   16: Solve the linear system

{I+δt(A−

Af ^((m)) , f ^((m))

)}f ^((m+1)) =f ^((m)).

-   -   17: Compute the mass sphere center c in (49);     -   18: Normalize f^((m+1)) by (f^((m+1))−c)/f^((m+1))−c∥ such that         the mass center of the resulting f^((m+1)) is at the unit sphere         center;     -   19: Compute δt according to strategies (S₂.a), (S₂.b) and (S₂.c)         with α=1.05 and ε=0.05;     -   20: until |ε_(h)(f^((m+1)))−ε_(h)(f^((m)))|<δΕ

Numerical Results

The performances of the computation of Riemann and sphere conformal mappings are evaluated by numerical experiments. The benchmark problems come from the eight different facial expressions for a person in FIG. 1 and four different genus zero closed surfaces in FIG. 2, respectively. The corresponding mesh data are listed in Table 1.

All computations are carried out in MATLAB 2011b on a HP workstation with two Intel Quad-Core Xeon X5687 3.6 GHz CPUs, 48 GB main memory and RedHat Linux operation system, using IEEE double-precision floating-point arithmetic.

1. Efficiency of Computing Riemann Mapping by Algorithm 1

The efficiency of Algorithm 1 in solving the nonlinear heat diffusion equation (10) with zero mass center by (i) comparing the efficiency of the explicit, implicit and quasi-implicit Euler methods, (ii) robustness in choosing initial time step and (iii) high performance computing for the benchmark human faces in FIG. 1, is demonstrated. In (i) and (ii), the human face in FIG. 1(a) is used as the benchmark problem. The stopping tolerance δΕ is taken as 10⁻⁹. The initial map f⁽⁰⁾ are all constructed by the harmonic-type initial map (see Section 5.1). Such f⁽⁰⁾ is close to the steady state solution so that the phase-I QIEM, i.e., lines 7-13, is skipped in Algorithm 1 Only phase-II QIEM, i.e., lines 15-20, is used to compute the steady state solution.

-   Comparison for the explicit, implicit and quasi-implicit Euler     methods: -   The time steps δt for the explicit, implicit and quasi-implicit     Euler methods are set to 0.003, 5 and 2881.05, respectively. The     energy ε_(h)(f^((m)) produced by Algorithm 1 is strictly     monotonically decreasing so that δt is fixed at each iteration. In     order to accelerate the convergence of the implicit Euler method, an     adaptive method is applied to control time step. For explicit Euler     method, δt is fixed in each iteration. The numerical results of the     explicit, implicit and quasi-implicit Euler methods are shown in     Table 2. The notation #it in Table 2 denotes the total iteration     number used to compute the solution. The convergence behaviors for     these three methods are shown in FIG. 3. From these numerical     results, it is learned that the QIEM outperforms the explicit and     implicit Euler methods greatly.

TABLE 1 #V and #F denote the numbers of vertices and faces of M, respectively, dim is the dimension of the matrix A. FIG. #F #V dim 1(a) 102256 51596 102258 1(b) 135766 68430 135768 1(c) 136513 68802 136515 1(d) 127265 64156 127267 1(e) 127132 64095 127134 1(f) 130560 65802 130562 1(g) 140352 70720 140354 1(h) 139198 70180 139200 2(a) 1005146 502575 502575 2(b) 858180 429092 429092 2(c) 284692 142348 142348 2(d) 39698 19851 19851

TABLE 2 Efficiency comparison Explicit Euler Implicit Euler QIEM δt 0.003 5 2,881.05 #it 2,412,000 108 10 CPU time (sec.) 113,911.89 5,554.88 4.23

-   Robustness of choosing initial time step in Algorithm 1: The     numerical results have validated that the QIEM with the time step     δt=2881.05, which is estimated by Algorithm 1, outperforms the other     two methods by big margins. Now, the robustness of Algorithml with     the harmonic-type initial map f (0) for solving (10) with human face     in FIG. 1(a) is demonstrated. Four different initial time steps δt₀,     namely (1, 10, 100, 1000), are applied to solve the steady state     problem and the maximum δt_(max) of time steps is equal to 2881.05.     The associated numerical results are shown in Table 3 and the     convergence behaviors of the harmonic energy are shown in FIG. 4.     The notation δt_(end) in Table 3 denotes the time step at the final     iteration. The results in Table 3 and FIG. 4 show that no matter     which initial time step is used, the adaptive controlling processes     in Algorithm 1 performs well. The time step in each iteration is     effectively increasing and the convergence of the solution is     accelerated. Furthermore, these results also show that the QIEM with     a harmonic-type initial f⁽⁰⁾ has a wide stability region and the     process of estimating time step in Algorithm 1 provides a near     optimal initial time step.

The adaptive controlling time step is not only robust for the human face problems but also for other benchmark problems. FIG. 5(b) presents the numerical results for applying Algorithm 1 to compute the Riemann mapping of the human hand in FIG. 5(a). As shown in FIG. 5(b), the time step is reduced from 3348.76 to 837.19 at the fourth iteration and then gradually increased to 3348.7. At the 41th iteration, it is again reduced to 1518.71 so that the harmonic energies are always decreasing.

-   In each iteration, the time step is kept greater than or equal to     837.19 such that it only requires 56 iterations and 103.3 seconds to     compute the Riemann mapping.

TABLE 3 Numerical results for Algorithm 1 with different initial time step δt₀. δt₀ 1 10 100 1000 δt _(end) 467.5 495.6 579.2 1551 #it 129 83 39 13 CPU time 52.4 33.9 16.1 5.63 (sec.)

-   High performance computing for Algorithm 1: the high performance     computing for Algorithm 1 with harmonic-type initial map f⁽⁰⁾ in     solving the benchmark problems in FIG. 1 is demonstrated. The     numerical results are shown in Table 4. The notation t_(c) in the     table denotes the total CPU time for computing the conformal map by     Algorithm 1 From Table 4, Algorithm 1 provides a large initial time     step δt₀. Using such δt₀, it only needs 7˜15 iterations in Algorithm     1 to compute the steady state solution of (10). The associated CPU     times are less than 17 seconds for FIG. 1, which illustrates the     high performance of Algorithm 1.

TABLE 4 Numerical results for the benchmark problems produced by Algorithm 1 with harmonic-type initial map f (0). FIG. δt₀ #it t_(c) (sec.) 1(a) 2881.1 10 4.23 1(b) 3377.7 8 4.80 1(c) 3375.1 7 4.34 1(d) 2890.6 11 6.02 1(e) 2155.6 10 5.46 1(f) 3307.4 15 8.11 1(g) 2102.3 13 8.08 1(h) 3474.0 11 16.9

2. Efficiency of Algorithm 1 for Computing Spherical Conformal Mapping

The high performance of Algorithm 1 in computing Riemann mapping has been numerically illustrated. In this subsection, the importance of phase-I QIEM in computing spherical conformal mapping is demonstrated and then the performance for Algorithm 1 in solving the benchmark problems in FIG. 2 is shown.

-   Effect of phase-I QIEM: Compare the efficiency of the phase-II QIEM     with Algorithm 1 for solving (10) with genus zero closed surface in     FIG. 2(c). The numerical results produced by Algorithm 1 are shown     in FIG. 6. Since the initial map f⁽⁰⁾ is far away from the steady     state solution, in numerical experiences, it is difficult to find an     initial time step such that the phase-II QIEM can be convergent. The     stability region of the time step in the phase-II QIEM for computing     spherical conformal mapping is extremely limited. On the contrary,     Algorithm 1 with initial time step δt₀.=2077.71 only requires 41     iterations which obviously outperforms phase-II QIEM. That is the     phase-I QIEM in Algorithm 1 plays an important role in computing     spherical conformal mapping.

TABLE 5 Numerical results for the benchmark problems in FIG. 2 produced by Algorithm 1 FIG. #it_(I) #it_(II) t_(c) (sec.) 2(a) 16 42 145.7 2(b) 8 121 234.6 2(c) 8 33 31.03 2(d) 5 80 6.117

-   Performance of Algorithm 1: the high performance computing for     Algorithm 1 in computing spherical conformal mapping for the     benchmark problems in FIG. 2 is demonstrated. The numerical results     are shown in Table 5. The notations #it_(I) and #it_(II) denote the     total iteration numbers for phase-I and phase-II QIEMs,     respectively. From the results in this table, it is concluded that     Algorithm 1 possesses high performance for computing the spherical     conformal maps of genus zero closed surfaces.

In summary, a mapping is conformal if and only if it is harmonic for genus zero closed surfaces. A traditional way to find the harmonic map is to minimize the harmonic energy by the time evolution according to the nonlinear heat diffusion (10). For this, an efficient quasi-implicit Euler method (QIEM) is proposed and applied it to compute conformal mappings from a genus zero closed surface to the unit sphere

² (the spherical conformal mapping) as well as from a simply connected surface with a single boundary to a 2D disk D (the Riemann mapping). Furthermore, the convergence of the QIEM under some simplifications is analyzed. In order to accelerate the convergence of this method, a variant time step scheme and a heuristic method to determine an initial time step have been developed. Numerical results validate that the proposed algorithms possess high performance for computing the Riemann mapping. For the spherical conformal map, a two-phase QIEM is proposed. In phase I QIEM, the normal component of Δ_(d)f is omitted to accelerate the computed map close to the steady state solution as soon as possible. When the computed map is close to the steady state solution in phase-I QIEM, it is switched to the adaptive time step controlling phase-II QIEM. Numerical results confirm that the two-phase QIEM also possesses high performance for computing the spherical conformal map.

Additional advantages and modifications will readily occur to those skilled in the art. Therefore, the invention in its broader aspects is not limited to the specific details, and representative devices shown and described herein. Accordingly, various modifications may be made without departing from the spirit or scope of the general inventive concept as defined by the appended claims and their equivalents. 

What is claimed is:
 1. A method for computing spherical conformal and Riemann mappings comprising the steps of: (a) carrying out evolution of computing a conthnnal map f from a genus zero closed surface to a unit sphere as well as from a simply connected surface with a single boundary to a 2D disk by a nonlinear heat diffusion equation; (b) solving the nonlinear heat diffusion equation by a quasi-implicit Euler method (QIEM); (c) analyzing convergence of the QIEM under some simplifications; (d) accelerating the convergence of the QIEM by using a two-phase approach for the quasi-implicit Euler method to estimate an initial time step and an adaptive method to control the time step.
 2. The method as claimed in claim 1, wherein a way to find the conformal map f is by time evolution according to the nonlinear heat diffusion equation in the step (a) ${\frac{df}{dt} = {{- \Delta_{d}}{f.}}};$ the time evolution of the conformal map f is according to the nonlinear heat diffusion equation $\frac{{df}(v)}{dt} = {{- \left( {\Delta_{d}{f(v)}} \right)^{}} = {- {\left( {{{{\Delta_{d}{f(v)}} -} < {\Delta_{d}{f(v)}}},{{n\left( {f(v)} \right)} > {n\left( {f(v)} \right)}}} \right).}}}$
 3. The method as claimed in claim 1, wherein the Quasi-Implicit Euler Method (QIEM) in the step (b) is iteratively formulated by $\begin{matrix} {\frac{f^{({m + 1})} - f^{(m)}}{\delta \; t} = {- \left( {{{{\Delta_{d}f^{({m + 1})}} -} \prec {\Delta_{d}f^{(m)}}},{f^{(m)} \succ f^{({m + 1})}}} \right)}} \\ {{= {{- \left( {{{A -} \prec {A\; f^{(m)}}},{f^{(m)} \succ}} \right)}f^{({m + 1})}}},} \end{matrix}$ wherein an unknown vector f^((m+1)) of F(f^((m+1)), f^((m)) ) is solved by Newton's method vec(f _(i+1) ^((m+1)))=vec(f _(i) ^((m+1)))−J(f _(i) ^((m+1)))⁻¹ F(f _(i) ^((m+1)) , f ^((m))) with f₀ ^((m+1))=f^((m)), wherein J(f^((m+1))) is the Jacobian matrix of F(f^((m+1)), f^((m))); wherein a new vector f^((m+1)) is generated in each iteration by solving a linear system [I+δt(A−

Af ^((m)) , f ^((m))

)]f ^((m+1)) =f ^((m)).
 4. The method as claimed in claim 1, wherein the two-phase approach for the quasi-implicit Euler method in the step (d) includes a phase-I QIEM formulated by ${\frac{{f^{({m + 1})} - f^{(m)}}}{\delta \; t} = {{{- \Delta_{d}}f^{({m + 1})}} = {{- A}\; f^{({m + 1})}}}},$ if f⁽⁰⁾ is not close to the steady state solution; i.e., (I+(δt)A)f ^((m+1)) =f ^((m)) is used to compute f^((m+1)); and a phase-II QIEM when f^((m)) is close to the steady state solution, switch to [I+δt(A−

Af^((m)), f^((m))

)]f^((m+1))=f^((m)). with repetitive strategies (S₂.a), (S₂.b) and (S₂.c) fir computing f^((m+1)) until difference |ε_(h)(f^((m+1)))−ε_(h)(f^((m)))| of energy is less than a given tolerance ε₃; wherein (S₁.a) By definition of A in ${a_{ii}:={\sum\limits_{{\lbrack{v_{i},v_{j}}\rbrack} \in M}\; k_{ij}}},{a_{ij}:={- k_{ij}}},{i \neq j},,$ row sums of A are equal to zero which implies that there is a trivial solution f such that (A−

Af, f

)f=0 and ε_(h)(f)=0; to avoid f^((m)) convergent to this trivial solution, reset δt₀ to be the current δt and restart algorithm with the original f⁽⁰⁾ and new δt₀ when ε_(h)(f^((m))) is less than a given small tolerance ε₂; (S₁.b) If ε_(h)(f^((m+1))) does not decrease, then δt is reduced to be δt:=max(1, 1/2 δt) and recompute f^((m+1)) from the [I+δt(A−

Af^((m)), f^((m))

)]f^((m+1))=f^((m)); (S₁.c) If the difference of the energy is still larger than or equal to the given tolerance ε₃ when m has been larger than a given maximal iteration number m_(max), reset δt₀ to be the current δt and restart the algorithm with original initial closed surface f⁽⁰⁾ and new δt₀. 